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Abstract 

In this article I review two data analysis methods for determining the mass (and even- 
tually the spin-independent cross section on nucleons) of Weakly Interacting Massive Par- 
ticles with positive signals from direct Dark Matter detection experiments: a maximum 
likelihood analysis with only one experiment and a model-independent method requiring 
at least two experiments. Uncertainties and caveats of these methods will also be discussed. 



1 Introduction 



There is strong evidence that more than 80% of all matter in the Universe is dark (i.e., interacts 
at most very weakly with electromagnetic radiation and ordinary matter) . The dominant compo- 
nent of this cosmological Dark Matter should be due to some yet to be discovered, non-baryonic 
particles. Weakly Interacting Massive Particles (WIMPs) x arising in several extensions of the 
Standard Model of electroweak interactions are one of the leading candidates for Dark Matter. 
WIMPs are stable particles with masses roughly between 10 GeV and a few TeV and interact 
with ordinary matter only weakly (for reviews of WIMPs and some other possible candidates 
for Dark Matter, see Refs. [1, 2, 3]). 

Currently, the most promising method to detect different WIMP candidates is the direct 
detection of the recoil energy deposited in a low-background laboratory detector by elastic 
scattering of ambient WIMPs on the target nuclei [4, 5, 6]^. The basic expression for the 
differential event rate for elastic WIMP-nucleus scattering is given by [1]: 
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Here R is the direct detection event rate, i.e., the number of events per unit time and unit mass 
of detector material, Q is the energy deposited in the detector, F{Q) is the elastic nuclear form 
factor, fi{v) is the one-dimensional velocity distribution function of the WIMPs impinging on 
the detector, v is the absolute value of the WIMP velocity in the laboratory frame. The constant 
coefficient A is defined as 

where po is the WIMP density near the Earth and (Tq is the total cross section ignoring the form 
factor suppression. The reduced mass ni^ n is defined by 



A^-^^, (2) 



mr,N = ^ , 3) 

where is the WIMP mass and mN that of the target nucleus. Finally, v^nin is the minimal 
incoming velocity of incident WIMPs that can deposit the energy Q in the detector: 



Vmin = ayQ (4) 

with 



and Umax is related to the escape velocity from our Galaxy at the position of the Solar system, 

^esc- 

It was found that, by using a time-averaged recoil spectrum dR/dQ, and assuming that no 
directional information exists, the normalized one-dimensional velocity distribution function of 
incident WIMPs, fi{v), can be solved from Eq.(l) directly as [7] 




(6) 



^Remind that, besides many different candidates for WIMPs, it is also possible that some other particles are 

(theoretically) candidates for Dark Matter. For more details about these various possible Dark Matter particles 
in many different (exotic) models or scenarios as well as the possible methods to detect them, see e.g., articles 
in Parts 1, 2, and 4 of this focus issue. 



where the normahzation constant H is given by 
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dQ . (7) 



Note that, firstly, because fiiy) in Eq.(6) is the normalized velocity distribution, the normahza- 
tion constant J\f here is independent of the constant coefficient A defined in Eq.(2). Secondly, 
the integral in Eq.(7) goes over the entire physically allowed range of recoil energies: starting at 
Q — 0, and the upper limit of the integral has been written as oo. However, it is usually assumed 
that the WIMP flux on the Earth is negligible at velocities exceeding the escape velocity v^sc- 
This leads thus to a kinematic maximum of the recoil energy 

QmaXjkin • (8) 

The velocity distribution function of halo WIMPs reconstructed by Eq.(6) is independent 
of the local WIMP density po well as of the WIMP-nucleus cross section (Tq. However, not 
only the overall normalization constant A/" given in Eq.(7), but also the shape of the velocity 
distribution, through the transformation Q = ja^ in Eq.(6), depends on the WIMP mass 
involved in the coefficient a defined in Eq.(5). In fact, any (assumed) value of will lead to 
a well-defined, normalized distribution function j\{v) when one uses Eq.(6). Hence, can be 
extracted from a single recoil spectrum only if one makes some assumptions about the velocity 
distribution fi{v). In contrast, by comparing two (or more) velocity distributions reconstructed 
from different recoil spectra with different target nuclei, one could avoid using these assumptions 
and estimate the WIMP mass model-independently. 

The remainder of this article is organized as follows. In Sec. 2 I first review a method for 
determining the WIMP mass with only one direct detection experiment. In Sec. 3 I present 
a model-independent method for determining by combining two experimental data sets. 
Numerical results based on Monte Carlo simulations of future experiments and uncertainties 
and caveats of these two methods will also be discussed. I conclude in Sec. 4. Some technical 
details for the data analysis will be given in an appendix. 



2 With one experiment 

In this section I review the method for determining the WIMP mass with only one direct detec- 
tion experiment based on a maximum likelihood analysis [8, 9, 10, 11]. 

2.1 Maximum likelihood analysis 

I first describe briefly some (standard) theoretical models/assumptions for fltting the elastic 
WIMP-nucleus scattering spectrum to experimental data. Then I discuss the determination of 
the WIMP mass by a maximum likehhood analysis. Note here that only the most commonly used 
models/assumptions are described as examples to show which information is required for the 
maximum likelihood analysis; however, it should be understood that other models or assumptions 
can also be used. 



2.1.1 Simple model distributions 



The simplest semi-realistic model halo is a Maxwellian halo. The one-dimensional velocity 
distribution function in the rest frame of our Galaxy can be expressed as [6, 1, 7] 
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Here Vq ~ 220 km/s is the orbital velocity of the Sun in the Galactic frame, and 
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is the normalization constant which satisfies 



Mv)dv^l. 
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Note that the second term on the right-hand side of Eq.(9) has been introduced to keep the 
velocity distribution continuous at v = v^sc- Substituting Eq.(9) into Eq.(l), the integral over 
the velocity distribution function can be calculated as 



(12) 

oo limit, iVcau ^ 4/\/^^o 



where v^an = OL\fQ in Eq.(4) has been used. Note that, in the fesc 
and the integral approaches to (2/y^fo) ^1'"^. 

On the other hand, when we take into account the orbital motion of the Solar system around 
the Galaxy as well as that of the Earth around the Sun, the velocity distribution function should 
be modified to [6, 1, 7] 
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for V < Vesc, with the normalization constant 
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Here Ve is the Earth's velocity in the Galactic frame [5, 1, 2]: 
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tp ~ June 2nd is the date on which the velocity of the Earth relative to the WIMP halo is 
maximal. Consequently, an analytic form of the integral over this velocity distribution can be 
given as 

rVesc r f , , 
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For practical, numerical uses, an approximate form of the integral over fi{v) was introduced as 
[6] 



where Cq and Ci are two fitting parameters of order unity. Not surprisingly, their values depend 
on the Galactic orbital and escape velocities, the target nucleus, the threshold energy of the 
experiment, as well as on the mass of incident WIMPs. Note that, the characteristic energy 

Qch ^ '-4 (18) 

and thus the shape of the recoil spectrum depend highly on the WIMP mass: for light WIMPs 
{m^ <^ rriTSi), Qch oc and the recoil spectrum drops sharply with increasing recoil energy, 
while for heavy WIMPs (m^ » rri^), Qch ~ const, and the spectrum becomes flatter. 

2.1.2 Local WIMP density 

Currently, the most commonly used value for the local WIMPs density in Eq.(2) is given as [1, 2] 

po ~ 0.3 GeV/cm^ . (19) 

However, so far it can be estimated only by means of the measurement of the rotational velocity 
of our Galaxy. Due to our location inside the Milky Way, it is more difficult to measure the 
accurate rotation curve of our own Galaxy than those of other galaxies. Thus an uncertainty of 
around a factor of 2 has been usually adopted [1, 2]: 

Po = 0.2 - 0.8 GeV/cm^ . (20) 

2.1.3 Spin-independent WIMP-nucleus cross section 

In most theoretical models, the spin-independent (SI) WIMP interaction on a nucleus with an 
atomic mass number ^4^30 dominates the spin-dependent (SD) interaction [1, 2]. Additionally, 

for the lightest supersymmetric neutralino, which is perhaps the best motivated WIMP candidate 
[1,2], and for all WIMPs which interact primarily through Higgs exchange, the SI scalar coupling 
is approximately the same on both protons p and neutrons n. The "pointlike" cross section ctq 
in Eq. (2) can thus be written as 
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and /p is the effective XXPP four-point coupling, A is the atomic mass number of the target 
nucleus. 



2.1.4 Nucleeir form factor 



For the SI cross section, an analytic nuclear form factor can be used. The simplest one is the 
exponential form factor, first introduced by Ahlen et al. [12] and Preese et al. [5]: 

F^iQ) = . (23) 

Here Q is the recoil energy transferred from the incident WIMP to the target nucleus, 
1.5 



is the nuclear coherence energy and 
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(24) 



(25) 



is the radius of the nucleus. The exponential form factor implies a Gaussian form of the radial 
density profile of the nucleus. This Gaussian density profile is simple, but not very realistic. 
Engel has therefore suggested a more accurate form factor [13], inspired by the Woods-Saxon 
nuclear density profile [1, 2], 



3ji(gi?i 
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Here ji{x) is a spherical Bessel function, 

q = \j2m^Q 
is the transferred 3-momentum, 



(26) 



(27) 



Ri = /R^^r5^ (28) 
is the effective nuclear radius^ with ^ 

^ 1.2A^/^ fm, (32) 

and 

s ~ 1 fm (33) 

is the nuclear skin thickness. 

^In the literature, the form factor given in Eq.(26) is also known as the "Helm" form factor with [14, 6] 

R, = ^ Rl+ (1)^2,2 -5S\ (29) 

where 

i?A (1.23 A^/^ -0.6) fm, ro =i 0.52 fm, s ~ 0.9 fm. (30) 

^For i?i given by Eq.(28) with s ~ 1 fm, a more precise approximation for Ra has also been given [15, 6]: 
Ra ^ (1.15 A^/^ + 0.39) fm. (31) 



2.1.5 Extended likelihood function 



Now we are ready to put all pieces for predicting the elastic WIMP-nucleus scattering spectrum 
together and then fit this spectrum to experimental data by maximizing the logarithm of the 
extended likehhood function [10]: 

Here 

\ = 8 \^\dQ (35) 



is the expected event number with the (assumed) exposure of the experiment, A^tot is the 
total number of events recorded in one (simulated) experiment, Qa are measured recoil energies 
in the data set between the minimal and maximal cut-off energies, Qmin and Qmax, and 

|-Qmax I dR\ 

is the total event rate. 

Note that, firstly, the definition of C in Eq.(34) takes into account the fact that the event 
number A'tot and the measured recoil spectrum S{dR/dQ) of each (simulated) experiment are 
not fixed. Secondly, except cq and ci in Eq.(17), there are two fitting parameters in the extended 
likelihood function £, i.e., the WIMP mass m-^ (involved in a) and the SI WIMP-proton cross 
section a^. 



2.2 Numerical results 

Here I show some numerical results with 10,000 simulated experiments based on Monte Carlo 
simulations performed by A. Green [10, 11]. ^^Ge has been chosen as the target nucleus with a 
threshold energy of 10 keV. A three-dimensional Maxwellian velocity distribution in the Galactic 
rest frame for an isotropic isothermal WIMP halo, taking into account the Earth's motion around 
the Sun with Vq = 220 km/s and fcsc = 540 km/s, and the Helm form factor in Eqs.(26), (27), 
(29), and (30) have been used. The standard assumption for the local WIMP density of 0.3 
GeV/cm^ has been adopted. 

Note that the simulations demonstrated here as well as in the next section for the method 
combining two experimental data sets are based on several simplified assumptions^. Firstly, the 
sample to be analyzed contains only signal events, i.e., is free of background. Active background 
suppression techniques [16, 17, 18]^ should make this condition possible. Secondly, all exper- 
imental systematic uncertainties as well as the uncertainty on the measurement of the recoil 
energy have been ignored. The energy resolution of most existing detectors is so good that its 
error can be neglected compared to the statistical uncertainty for the foreseeable future. 



2.2.1 Statistical uncertainty 

Figs. 1 show the distributions of the best-fit WIMP mass and SI WIMP-proton cross section 
(T^p on the cross section versus WIMP mass plane. The input WIMP mass and the cross 

^More realistic modelling with e.g., other WIMP velocity distributions and/or different nuclear form factors 

could in principle be incorporated into the maximum likelihood analysis. 

^ For more experimental details about current direct detection techniques and the next generation detectors, 
see articles in Part 3 of this focus issue. 
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Figure 1: Distributions of the best-fit WIMP mass and SI WIMP-proton cross section on the 
cross section versus WIMP mass plane. The input WIMP mass and the cross section are 100 
GeV and 10^'' pb, respectively. The exposures have been assumed to be 3 x 10^ (left) and 3 x 10^ 
(right) kg-day and the corresponding expected event numbers are 78 and 780, respectively. In 
each frame, the contours contain 68% and 95% of the simulated experiments. See the text for 
further details (Plots from [10]). 



section are 100 GeV and 10^^ pb, respectively. The exposures have been assumed to be 3 x 10^ 
(left) and 3 x 10^ (right) kg-day and the corresponding expected event numbers are 78 and 
780^, respectively. It can be seen that, especially for the smaller exposure, the distribution 
is asymmetric and there are (significantly) more experiments with best-fit masses and cross 
sections larger than the inpiit values. QTiantitativcly, for a WIMP mass of 100 GeV with ~ 80 
events, the la and 2a statistical uncertainties are I35 GeV and GeV, respectively [10]. 

Fig. 2 shows the 95% (solid) and 68% (dotted) confidence limits on the best-fit WIMP mass 
as functions of the input WIMP mass. The input SI WIMP-proton cross section has been set 
here as 10~^ pb. The assumed exposures are 3 x 10^, 3 x 10^, and 3 x 10^ kg-day, respectively. 
We see here that since, as mentioned above, the shape of the recoil spectrum varies significantly 
with the WIMP mass for light WIMP masses (m^ < m^), the WIMP mass (and also the cross 
section) can be fitted with a higher accuracy: the la and 2a statistical uncertainties for = 25 
GeV are ±4 GeV and GeV, for = 50 GeV are t\l GeV and GeV, respectively [10]. 

In contrast, the weak dependence of the shape of the recoil spectrum on the WIMP mass 
for heavy WIMP masses (m^ ^ ttin) means that it will be more difficult or even impossible 
to extract the WIMP mass with (9(100) events, if WIMPs are (much) heavier than the target 
nucleus [10]. Note that the dependence of the shape of the recoil spectrum on the WIMP mass 
as well as on that of the target nucleus suggests that heavy nuclei, e.g., Xe, would be able to 
measure the mass of heavy WIMPs more accurately; however, the rapid decrease of the nuclear 
form factor with increasing recoil energy, which occurs for heavy nuclei, means that, due to less 
expected ewnts. this is in fact not necessarily the case. 



^ Since the event number is directly proportional to the product of the cross section and the exposure £ 



it is equivalent to assume o-yi, = 10 ° pb and exposures of 3 x 10 and 3 x 10° kg-day. 




50 100 150 200 

(GeV) 

Figure 2: The 95% (solid) and 68% (dotted) confidence limits on the best-fit WIMP mass as 
functions of the input WIMP mass. The input SI WlMP-proton cross section has been set here 
as 10~^ pb. The assumed exposures are 3 x 10^, 3 x 10^, and 3 x 10^ kg-day, respectively (Plot 
from [11]). 
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Figure 3: Distributions of the best-fit WIMP mass and SI WlMP-proton cross section on the 
cross section versus WIMP mass plane. The input orbital velocity of the Solar system vo has 
been set as 200 (left) and 240 (right) km/s, while the standard value of vo — 220 km/s has been 
used for the data analysis. The exposure assumed here is 3 x 10^ kg-day. The other parameters 
are as in Figs. 1 (Plots from [10]). 



2.2.2 Systematic uncertainties 



Different sources of the systematic uncertainties in this model-dependent analysis have been 
considered [10, 11]. Figs. 3 show the distributions of the best-fit WIMP mass and cross section 
with different input orbital velocity of the Solar system: Vo = 200 (left) and 240 (right) km/s, 
while the standard value of Vq = 220 km/s has been used for the data analysis. As shown here, 
for an input WIMP mass of 100 GeV, there could be an ±20 GeV shift in the best-fit WIMP 
mass combined with an ~ ±10"*^ pb (~ 10%) shift in the SI WIMP-proton cross section caused 
by the ± 20 km/s difference between the real and the assumed orbital velocities [10]. Moreover, 
the larger the real orbital velocity, the less the expected event number (with a fixed exposure), 
and thus the larger the statistical uncertainties on both the WIMP mass and SI WIMP-proton 
cross section one could obtain. 

More detailed illustrations and discussions about the effects of varying the underlying WIMP 
mass and cross section, the detector target nucleus, the exposure, the minimal and maximal cut- 
off energies, the orbital velocity of the Solar system, as well as the background event rate and 
its spectrum can be found in Refs. [10, 11, 19]. 



3 Combining two experiments 

In this section I first review the model-independent method for reconstructing the WIMP mass 
by using two experimental data sets with different target nuclei^. Then I also describe an 
extension of this method for estimating (or at least constraining) the SI WIMP-proton cross 
section. 



3.1 Model— independent determination 

As mentioned in the introduction, the normalized one-dimensional velocity distribution function 
of incident WlMPs can be solved from Eq.(l) directly and, consequently, its generalized moments 
can be estimated by [20] 

rv(Qina;x) 

(t;")(t;(gmin),w(Qmax)) = / ,„ , v"" fi{v) dv 

(37) 
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Here v{Q) = ay/Q, (5(min,max) are the experimental minimal and maximal cut-off energies, 

r(Q.in) ^ (^] (38) 

V'^^/ expt, Q=Q^in 

is an estimated value of the measured recoil spectrum {dR/dQ)cxpt {before the normalization by 
the exposure S) at Q — Qmm, and In{Qmm, Qmax) can be estimated through the sum: 

g(n-l)/2 



''In Ref. [8], the authors mentioned an attempt for using the maximum Ukelihood analysis with two (or more) 

detector materials. However, they found that, since the likelihood contours for different targets are pretty similar 
when simulating with the same number of events, their results showed effectively little different from that obtained 
with a single experiment. 



where the sum runs over all events in the data set that satisfy Qa £ [Qmin, <5max]- Note that, 
firstly, by using the second Eq.(37) (v")(t'((5min), t'((5niax)) can be determined independently of 
the local WIMP density po) of the velocity distribution function of incident WlMPs, ), as well 
as of the WlMP-nucleus cross section ctq. Secondly, as shown later, r((5mm) and IniQmim Qmax) 
are two key quantities for this model-independent method, which can be estimated either from 
a functional form of the recoil spectrum or from experimental data (i.e., the measured recoil 
energies) directly^. However, r((5min) and /n(Qminj Qmax) estimated from a scattering spectrum 
fitted to experimental data are not model-independent any more. 



3.1.1 Basic expressions for determining m. 



By requiring that the values of a given moment of fi{v) estimated by Eq.(37) from two detectors 
with different target nuclei, X and F, agree, my_ appearing in the prefactor a" on the right-hand 
side of Eq.(37) can be solved as [21]: 

_ ^mxiriY - mx{Tln,x/T^n,Y) 



m 



where 
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(40) 



'^Qmin,xfx(Qmin,x) / Fx(Qmin,x) + Io,X 



1/n 



(41) 



and TZnx can be defined analogously. Here n ^ 0, m(^x,Y) -^(x,y)(Q) are the masses and the 
form factors of the nucleus X and Y, respectively, and r(x,y)((5min,(x,y)) refer to the counting 
rates for detectors X and Y at the respective lowest recoil energies included in the analysis. Note 
that, firstly, the general expression (40) can be used either for spin-independent or for spin- 
dependent scattering, one only needs to choose different form factors under different assumptions. 
Secondly, the form factors in the estimate of In,x and In,Y using Eq.(39) are also different. 

On the other hand, by using the theoretical prediction that the SI WIMP-nucleus cross 
section given in Eq.(21) dominates, and the fact that the integral over the one-dimensional 
WIMP velocity distribution on the right-hand side of Eq.(l) is the minus-first moment of this 
distribution, which can be estimated by Eq.(37) with n~—l, one can easily find that [20] 
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Note that the exposure of the experiment, £^ appears in the denominator. Since the unknown 
factor Pol/pP on the left-hand side above is identical for different targets, it leads to a second 
expression for determining [20] 

[mx/myf^'^ my - mx{Tla,x lT^a,Y) 
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Here ?7^(x,y) has been assumed. 
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(44) 



and similarly for TZ^y. 



^AU formulae needed for estimating r((5min), I n{Q mm iQ max) ^ and their statistical errors are given in the 
appendix. 



3.1.2 x^-fi^tting 



In order to yield the best-fit WIMP mass as well as to minimize its statistical error by combining 
the estimators for different n in Eq.(40) with each other and with the estimator in Eq.(43), a 
function has been introduced [20] 
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the other ^max + 2 functions fi^y can be defined analogously. Here rimax determines the highest 
moment of fi{v) that is included in the fit. The fi are normalized such that they are dimension- 
less and very roughly of order unity in order to alleviate numerical problems associated with the 
inversion of their covariance matrix. Note that the first nmax + 1 fit functions depend on only 
through the overall factor a and that in Eqs. (46a) and (46b) is now a fit parameter, which 
may differ from the true value of the WIMP mass. Finally, C in Eq.(45) is the total covariance 
matrix. Since the X and Y quantities are statistically completely independent, C can be written 
as a sum of two terms^: 



dj = COV (Ax, fj,x) + COV (/j,y, Ay) . 



(47) 



3.1.3 Matching the cut— off energies 

The basic requirement of the expressions for determining given in Eqs. (40) and (43) is that, 
from two experiments with different target nuclei, the values of a given moment of the WIMP 
velocity distribution estimated by Eq.(37) should agree. This means that the upper cuts on 
fi{v) in two data sets should be (approximately) equaP°. Since v^ut — (^VQmax, it requires that 
[20] 



Q 



max,y 



ay/ 



(48) 



Note that a defined in Eq.(5) is a function of the true WIMP mass. Thus this relation for 
matching optimal cut-off energies can be used only if is already known. One possibility to 
overcome this problem is to fix the cut-off energy of the experiment with the heavier target, 
minimize the x^(m^) function defined in Eq.(45), and estimate the cut-off energy for the lighter 
nucleus by Eq.(48) algorithmically [20]. 



^Formulae needed for estimating the entries of C will be given in the appendix. 
^°Here the threshold energies have been assumed to be negligibly small. 



3.2 Numerical results 



Here I show some numerical results for the reconstructed WIMP mass based on Monte Carlo 
simulations. The upper and lower bounds on the reconstructed WIMP mass are estimated from 
the requirement that exceeds its minimum by l^-*^. ^^Si and ^^Gc have been chosen as two tar- 
get nuclei. The scattering cross section has been assumed to be dominated by spin-independent 
interactions. The shifted Maxwellian velocity distribution given in Eq.(13) (the second term in- 
volving Vesc has been neglected) with = 220 km/s, = l.Obvo^^, and v^sc ~ 700 km/s and the 
Woods-Saxon form factor in Eq.(26) have been used. The threshold energies of two experiments 
have been assumed to be negligible and the maximal experimental cut-off energies are set as 
100 keV. 2 X 5,000 experiments have been simulated. In order to avoid large contributions from 
very few events in the high energy range to the higher moments [7], only the moments up to 
iT-max — 2 were included in the fit. 



3.2.1 Statistical uncertainty 



In Figs. 4 the dotted (green) curves show the median reconstructed WIMP mass and its la- 
upper and lower bounds for the case that both Qmax.Si and (5max,Ge have been fixed to 100 keV. 
As argued earlier, the values of a given moment of the WIMP velocity distribution estimated by 
Eq.(37) do not agree when the same maximal cut-off energy for both experimental data sets is 
used. This causes a systematic underestimate of the reconstructed WIMP mass [21] which can 
be seen obviously here. 

The sohd (black) curves were obtained by using Eq.(48) for matching the cut-off energy 
Qmax,si perfectly with Qmax.Gc = 100 keV and the true (input) WIMP mass, whereas the dashed 
(red) curves show the case that Qmax.Gc = 100 keV, and Qmax.Si h^-s been determined by minimiz- 
ing x^(m^; Qmax.Si)- As shown here, with only 50 events on average before cuts (upper frame) 
from each experiment, the algorithmic process seems already to work pretty well for WIMP 
masses up to ~ 500 GeV. For ^ 100 GeV the median WIMP mass determined in this way 
overestimates its true vahie by 15 to 20%; however, the true WIMP mass always lies within the 
median limits of the la statistical error interval estimated by the algorithmic Qmax matching 
procedure up to even rriy = 1 TeV 



3.2.2 Statistical fluctuation 

In order to study the statistical fluctuation of the reconstructed WIMP mass by algorithmic 
Qmax matching in the simulated experiments, an estimator 5m has been introduced as [20] 
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(49) 



"'^"'^The median, rather than the mean, values for the (bounds on the) reconstructed WIMP mass are shown. 
^^The time dependence of the Earth's velocity in the Galactic frame, the second term of Ve{t) in Eq.(15), has 
been ignored. 
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Figure 4: Results for the reconstructed WIMP mass as well as its la statistical error interval 
based on the x^^fit in Eq.(45). 50 (upper) and 500 (lower) events on average before cuts from 
each experiment have been simulated. See the text for further details (Plots from Ref. [20]). 
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Figure 5: Normalized distribution of the estimator 6m defined in Eq.(49) for an input WIMP 
mass of 50 GeV with 50 events on average (before cuts) in each experiment. The other parameters 
and notations are as in Figs. 4 (Plot from Ref. [20]). 

Here m-^^in is the true (input) WIMP mass, rn-^^^ec its reconstructed value, m^joi{2) are the 1 (2) a 
lower bounds satisfying x^("^x,io(i,2)) = X^('^x,rec) + 1 (4), and m;^^hii(2) are the corresponding 
1 (2) a upper bounds. It has been found that the error intervals of the median reconstructed 
WIMP mass are quite asymmetric; similarly, the distance between the 2(T and la hmits can 
be quite different from the distance between the la limit and the central value [20]^^. The 
definition of 5m in Eq.(49) takes these differences into account, and also keeps track of the sign 
of the deviation: if the reconstructed WIMP mass is larger (smaller) than the true one, 5m is 
positive (negative). Moreover, \5m\ < 1 (2) if and only if the true WIMP mass lies between the 
experimental 1 (2) a limits. 

Fig. 5 shows the distribution of 5m calculated from 5,000 simulated experiments with 50 
events on average before cuts for a rather light WIMP mass of 50 GeV. In this case simply fixing 
both Qmax values to 100 keV still works fine (see the upper frame of Figs. 4). However, the 
distributions for both fixed Qmax and optimal Qmax matching look somewhat lopsided, since the 
error interval is already asymmetric, with TTi^ hii ~ "^x.rec > ""^x.^c ~ "^x.ioi- '^^^ overestimate of 
light WIMP masses reconstructed by algorithmic Qmax matching shown in Figs. 4 is reflected 
by the dashed (red) histogram here, which has significantly more entries at positive values than 
at negative values. These distributions also indicate that the statistical uncertainties estimated 
by minimizing ^(^{m^ are indeed overestimated, since nearly 90% of the simulated experiments 
have \5m\ < 1 [20], much more than ~ 68% of the experiments, that a usual la error interval 
should contain. 



Recall that the same asymmetry has also been observed by the maximum likelihood analysis. 
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Figure 6: Normalized distribution of the estimator Sm. Parameters and notations are as in 
Fig. 5, except that the input WIMP mass has been increased to 200 GeV. In the right frame the 
average event number (before cuts) in each experiment have, in addition, been increased from 50 
to 500. Note that the bins at 5m = ±5 are overflow bins, i.e., they also contain all experiments 
with \6m\>5. (Plots from Ref. [20]). 



Unfortunately, as shown in Figs. 6, when the true (input) WIMP mass increases to 200 GeV 
and the expected event number (before cuts) increases to 500 (right frame), the situations become 
less favorable. While optimal Qmax matching seems to approach very slowly to be Gaussian and 
the overestimated statistical errors become a little bit more reliable for larger event numbers [20], 
the errors estimated by the algorithmic procedure for determining Qmax.Si are not very reliable 
in the simulations. 

More detailed illustrations and discussions about algorithmic Qmax matching with different 
detector materials or with data sets generated in different halo models, as well as about the 
statistical fluctuation in the analysis can be found in Ref. [20]. 

3.3 Estimating the SI WIMP— proton coupling 

In the maximum likelihood analysis discussed in Sec. 2, the SI WIMP-proton cross section cj^p is 
the second fitting parameter that, combined with the WIMP mass m^, maximizes the extended 
likelihood function C calculated from an assumed WIMP velocity distribution. 

In contrast, as shown above, by combining two experimental data sets, one can estimate the 
WIMP mass without knowing the WIMP-nucleus cross section a^. Conversely, by means of 
Eq.(42), one can also estimate or at least constrain the SI WIMP-proton coupling, |/pp, from 
experimental data directly without knowing the WIMP mass [22]. 

3.3.1 Making an assumption for the local WIMP density 

In Eq.(42) the WIMP mass on the right-hand side can be determined by the method de- 
scribed above, r((5min) and Iq can also be estimated from one of the two data sets used for 
determining or from a third experiment. Nevertheless, due to the degeneracy between the 
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Figure 7: Left: The reconstructed SI WIMP-proton coupling l/pl^ec 
WIMP mass m 



as a function of the input 



Right: The reconstructed coupling |/p|^ec the reconstructed WIMP mass 



m 



on the cross section (coupling) versus WIMP mass plane. The open (red) squares indicate 
the input WIMP masses and the true values of the coupling. The open (blue) circles and the 
(blue) crosses indicate the reconstructed couplings and their la statistical errors. The horizontal 
and vertical solid (blue) lines show the la statistical errors on m^^rec and |/p|rec) respectively. 
Parameters are as in Figs. 4, in addition cr^p has been set as 10~^ pb. Each experiment has 50 
events on average. See the text for further details (Plots from Ref. [22]). 



local WIMP density po and the coupling |/pp, one cannot estimate both of them independently. 
The simplest way is making an assumption for the local WIMP density po^^- 



3.3.2 Numerical results 

The left frame of Figs. 7 shows the reconstructed SI WIMP-proton coupling |/p|rec as a function 
of the input WIMP mass m^^in. Following simulations for the reconstruction of the WIMP mass, 
^^Si and ''^Ge were chosen as two target nuclei for estimating in Eq.(42). In order to avoid 
complicated calculations of the correlation between the error on the reconstructed and that 
on the estimator of Jq, a second, independent data set with ™Ge was chosen as the third target 
for estimating Jq. The SI WIMP-proton cross section was set as 10~^ pb. Each experimental 
data set has 50 events on average under the common experimental cut-off energy Qmax chosen 
as 100 GeV. 

It can be seen that the reconstructed |/pp values are underestimated for WIMP masses ^ 100 
GeV. This systematic deviation is caused mainly by the underestimate of Jq. However, in spite 
of this systematic deviation (and in fact due to the fairly large statistical uncertainty), the true 
value of |/pp always lies within the la statistical error interval. Moreover, for a WIMP mass of 
100 GeV, one could in principle already estimate the SI WIMP-proton coupling with a statistical 
uncertainty of only ~ 15% with just 50 events from each experiment. Recall that this is much 
smaller than the systematic uncertainty of the local Dark Matter density (of a factor of 2 or 



^"'Note that, since the couphng |/pp estimated by Eq.(42) is inversely proportional to the local density po, 
whose common value falls on the lower end of the possible range (see Eqs.(19) and (20)), one can therefore at 
least give an upper bound on this coupling. 



even larger). 

Combining the estimate for the SI WIMP-proton couphng with the estimate for the WIMP 
mass, the right frame of Figs. 7 shows the reconstructed coupling |/p|rcc ^'^^ the reconstructed 
WIMP mass m^^rec on the cross section (coupling) versus WIMP mass plane^^. It is important to 
note that, as shown here, |/p|^ and can be estimated separately and from experimental data 
directly with neither prior knowledge of each other nor an assumption for the WIMP velocity 
distribution. 

4 Summary and conclusions 

In this article I reviewed the methods for the determination(s) of the mass (and eventually 
the spin-independent cross section on nucleons) of Weakly Interacting Massive Particles with 
positive signals of their elastic scattering off target nuclei in direct Dark Matter detection ex- 
periments. 

With only one experiment, the WIMP mass combined with its SI cross section on nucleons 
could be estimated by the maximum likelihood analysis using a theoretically predicted scattering 
spectrum fitted to the measured recoil energies. If WlMPs are light (m^ < m^), the shape of 
the recoil spectrum is sensitive to their mass, then the WIMP mass (and also the cross section) 
can be estimated with a higher accuracy; however, in case WIMPs are (much) heavier than the 
target nucleus (m^ ^ 200 GeV), the recoil spectrum becomes nearly independent on and it 
is then more difficult or even impossible to estimate the WIMP mass reasonably with C(IOO) 
events. 

The maximum likelihood analysis depends on the prior assumption for the velocity distri- 
bution of halo WIMPs as well as on the local WIMP density. For a WIMP mass of 100 GeV, 
an ~ 10% measurement uncertainty on the orbital velocity of the Solar system could cause an 
~ 20% systematic error on the best-fit WIMP mass combined with an ~ 10% error on the SI 
WIMP-proton cross section. 

In order to determine the WIMP mass without making any assumption for the WIMP ve- 
locity distribution, I described a second method based on the reconstruction of (the moments 
of) the WIMP velocity distribution function from two experiments with different target nuclei. 
This method can be used without knowing the WIMP-nucleus cross section. The only infor- 
mation needed is the measured recoil energies. By matching the maximal cut-off energies of 
two experiments one could in principle estimate the WIMP mass up to ~ 500 GeV with C(50) 
events from each experiment. 

Nevertheless, the algorithmic procedure for determining the maximal cut-off energy of the 
experiment with the lighter target nucleus by minimizing could overestimate the WIMP 
mass by 15 to 20% if WIMPs are light, or lead to unreliable error estimates if WIMPs are heavy. 
The latter could become worse with larger event samples. However, the fact that optimal 
Qmax matching works well in all cases, for both the median reconstructed WIMP mass and its 
statistical error, gives us hope that a better algorithm for Qmax matching can be found which 
only relies on the data. 

Additionally, by combining two (or three) experimental data sets one could also estimate the 
spin-independent WIMP-proton coupling without knowing the WIMP mass. Although, due to 
the degeneracy between the local WIMP density and the WIMP-nucleus cross section, one needs 
to adopt the local Dark Matter density (as the unique assumption), at least an upper bound on 



Plots shown here have been calculated by a different program than that for the Monte Carlo simulations 
shown in Figs. 4 to 6. 



this coupling could be given. In fact, for a WIMP mass of 100 GeV, with (9(50) events from 
each experiment, a statistical uncertainty of ~ 15% could be reached. This is much smaller than 
the systematic uncertainty on the local Dark Matter density (of a factor of 2 or even larger). 

In summary, by means of currently running and projected experiments using detectors with 
10~^ to 10~^^ pb sensitivities [16, 17, 18] (see footnote 5), we stand a good chance of detecting 
Dark Matter particles, if Dark Matter indeed consists (mainly) of WIMPs. Then the methods 
presented here can be used to estimate the mass (and eventually the cross section on nucleons) 
of Dark Matter particles. This information (perhaps combined with information from indirect 
detection experiments [19]) will allow us not only to constrain the parameter space in different 
extensions of the Standard Model of particle physics, but also to identify WIMPs among new 
particles produced at colliders (hopefully in the near future). Once one is confident of this 
identification, one can use further collider measurements of the mass and couplings of WIMPs. 
Together with the reconstruction of the velocity distribution of halo WIMPs [7], this will then 
yield a new determination of the local WIMP density. On the other hand, knowledge of the 
WIMP couplings will also permit prediction of the WIMP annihilation cross section. Together 
with information on the WIMP density, this will allow one to predict the event rate in the 
indirect Dark Matter detection [1, 2] as well as to test our understanding of the early Universe. 
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A Formulae needed in Sec. 3 

Here I list all formulae needed in the model-independent method described in Sec. 3. Detailed 
derivations and discussions can be found in Refs. [7, 20]. 

A.l Estimating r(Qrnin)5 -?^u(Qmin? Qmax), and their statistical errors 

First, consider experimental data described by 



Here the total energy range between Qmin and Qmax has been divided into B bins with central 
points Qn and widths bn- In each bin, Nn events will be recorded. Since the recoil spectrum 
dR/dQ is expected to be approximately exponential, the following ansatz for the spectrum in 
the nth bin has been introduced [7]: 
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Here 



r„ is the standard estimator for dR/dQ &X, Q — Qn- 



(A3) 



kn is the logarithmic slope of the recoil spectrum in the nth bin, which can be computed numer- 
ically from the average Q— value in the nth bin: 
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The error on the logarithmic slope kn can be computed from Eq. (A4) directly: 
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Qs^n in the ansatz (A2) is the shifted point at which the leading systematic error due to the 
ansatz is minimal [7], 
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Note that (5.s,„ differs from the central point of the nth bin, Qn- Prom the ansatz (A2), the 
counting rate at Q — Qmin can be calculated by 

r(Qn.in) = rie'=^(«-'"-«-i), (A9) 

and its statistical error can be expressed as 



<j'(r(gn.in)) = r2(Qmin) <^ TF + 



1 



— - — 1 + coth 



'biki 



<y\k. 



smce 



2/ \ 

'J (r„) = . 

Finally, since all are determined from the same data, they are correlated with 

Q{n+m-2)/2 



COv(/„, Im) = 



a F^Qa) ' 



(AlO) 



(All) 



(A12) 



where the sum again runs over all events with recoil energy between Qmin and Qmax- And the 
correlation between the errors on r((5inin), which is calculated entirely from the events in the 
first bin, and on is given by 
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note that the sums li here only count in the first bin, which ends at Q = Qmin + 

On the other hand, with a functional form of the recoil spectrum (e.g., fitted to experimental 

data), {dR/ dQ) cxpt, one can use the following integral forms to replace the summations given 
above. Firstly, the average Q— value in the nth bin defined in Eq.(A5) can be calculated by 
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For IniQram, Qmax) givcu in Eq.(39), we have 

rQmax g('^-i)/2 fdR\ 
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and similarly for the covariance matrix for /„ in Eq.(Al2), 
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Remind that {dR/dQ)epxt is the measured recoil spectrum before the normalization by the expo- 
sure. Finally, /^(Qmin, <5min + bi) needed in Eq.(Al3) can be calculated by 
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Note that r{Qrain) and In{Qmm,Qrain + &i) should be estimated by Eqs.(A9) and (A17) with 
ri, ki and Qs,i estimated by Eqs.(A3), (A4), and (A8) in order to use the other formulae for 
estimating the (correlations between the) statistical errors without any modification. 



A. 2 Statistical errors on given in Eqs.(40) and (43) 

The expression for m^\^^„^ given in Eq.(40) leads to a lengthy expression for its statistical error: 
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Here a short-hand notation for the six quantities on which the estimate of m-^ depends has been 
introduced: 



Ci,x — ^n,X , 



C2,X — Io,x , 



C3,X — '"x(Qmin,x) ; 



(A19) 



and similarly for the Qy. Estimators for cov(cj,Cj) have been given in Eqs.(Al2) and (A13). 
Explicit expressions for the derivatives of TZn,x with respect to Ci^x are: 
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explicit expressions for the derivatives of TZn,Y with respect to q y can be given analogously. Note 
that, firstly, factors TZn.(x.Y) appear in all these expressions, which can practically be cancelled 
by the prefactors in the bracket in Eq.(A18). Secondly, all the Io,x, -^o,y, In,x, In,Y should be 
understood to be computed according to Eqs.(39) or (A15) with integration limits Qmin and 
Qmax specific for that target. 

Similar to the analogy between Eqs.(40) and (43), the statistical error on my,\^ given in 
Eq.(43) can be expressed as 
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where we have again used the short-hand notation in Eq.(Al9); note that ci,(x,y) = In,{x,Y) 
not appear here. Expressions for the derivatives of TZa,x can be computed from Eq.(44) as 
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_'^Qmin,x''^x(Qmin,x) + -^0,X-^l(Qmin,x) 

and similarly for the derivatives of TZ^-x- 

A. 3 Covariance of fi defined in Eqs.(46a) and (46b) 

The entries of the C matrix in Eq.(47) involving basically only the moments of the WIMP 
velocity distribution can be read off Eq.(82) of Ref. [7], with an slight modification due to the 
normalization factor in Eq.(46a)^^: 

cov (/,, /,) = Af^ [/, cov(7o, lo) + Oi'^'{i + 1) (j + l)cov(7,, /,) 

- OL\i + cov(Jo, Ij) - a'{i + l)fj cov(/o, h) 

+ DiDja^riQ^^)) - {Difj + Djfi) cov{r{Q^^), Jq) 

+ a^{j + 1) A cov{r{Q^in), Ij) + a\i + l)Dj coy{r{Q^in)Ji) . 

(A23) 



^^Since the last fi defined in Eq.(46b) can be computed from the same basic quantities, i.e., the counting rates 
at Qmin and the integrals lo, it can directly be included in the covariance matrix. 
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A. 4 Statistical error on |/pp given in Eq.(42) 

From Eq.(42), it can easily be found that 
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where jVm is defined in Eq.(A24), and 
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The correlation between the error on the reconstructed and that on the estimator of l/TVm, 
the third term in Eq.(A27), can be neglected in case one uses three independent data sets. 
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